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The surface zonal winds observed in the giant planets form a complex jet pattern with alternating prograde and retrograde direction. 
While the main equatorial band is prograde on the gas giants, both ice giants have a pronounced retrograde equatorial jet. 

We use three-dimensional numerical models of compressible convection in rotating spherical shells to explore the properties of 
zonal flows in different regimes where either rotation or buoyancy dominates the force balance. We conduct a systematic parameter 
study to quantify the dependence of zonal flows on the background density stratification and the driving of convection. 

In our numerical models, we find that the direction of the equatorial zonal wind is controlled by the ratio of the global- scale 
buoyancy force and the Coriolis force. The prograde equatorial band maintained by Reynolds stresses is found in the rotation- 
dominated regime. In cases where buoyancy dominates Coriolis force, the angular momentum per unit mass is homogenised and 
the equatorial band is retrograde, reminiscent to those observed in the ice giants. In this regime, the amplitude of the zonal jets 
depend on the background density contrast with strongly stratified models producing stronger jets than comparable weakly stratified 
cases. Furthermore, our results can help to explain the transition between solar-like (i.e. prograde at the equator) and the "anti-solar" 
differential rotations (i.e. retrograde at the equator) found in anelastic models of stellar convection zones. 

In the strongly stratified cases, we find that the leading order force balance can significantly vary with depth. While the flow 
in the deep interior is dominated by rotation, buoyancy can indeed become larger than Coriolis force in a thin region close to the 
surface. This so-called "transitional regime" has a visible signature in the main equatorial jet which shows a pronounced dimple 
where flow amplitudes notably decay towards the equator. A similar dimple is observed on Jupiter, which suggests that convection 
in the planet interior could possibly operate in this regime. 

Keywords: Atmospheres dynamics, Jupiter interior, Saturn interior, Uranus interior, Neptune interior 
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1. Introduction 

Surface zonal jets on the giant planets Jupiter, Saturn, Uranus 
and Neptune have been investigated since the 70s by tracking 
cloud features. On each planet, these zonal winds form a differ- 
ential rotation profile with alternating prograde (i.e. eastward) 
and retrograde (i.e. westward) flows. 

On both gas giants, a strong eastward equatorial jet is flanked 
■by multiple weaker alternating zonal winds (10-20 m.s -1 ). 
Jupiter's equatorial jet extends roughly between ±20° latitude 
reaching a maximum velocity around 150 m.s -1 (Porco et al., 
2003; Vasavada and Showman, 2005). Saturn's equatorial jet is 
fiercer and wider with a maximum flow amplitude of 450 m.s" 1 
and extension of ±35° latitude (Porco et al., 2005). 

The zonal wind profiles of Uranus and Neptune are quite 
different. One broad retrograde equatorial jet is flanked by 
only two strong prograde jets at higher latitudes. On Uranus, 
this equatorial sub-rotation extends between ±30° latitude and 
reaches 100 m.s -1 (Hammel et al., 2005); it extends over ±50° 
latitudes and reaches 400 m.s -1 on Neptune (Sromovsky et al., 
1993). 
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Two competing categories of models try to explain the ob- 
served zonal winds structure. In the "shallow-forcing" sce- 
nario, zonal winds are driven by injection of turbulence via dif- 
ferent types of physical forcings at the stably- stratified cloud 
level (e.g., latent heat release, solar radiation, moist convec- 
tion). These models successfully reproduce an alternating zonal 
flow pattern similar to those observed in giant planets (e.g. 
Williams, 1978; Cho andPolvani, 1996). Although many ear- 
lier shallow models predict a similar westward equatorial zonal 
flow for all four giant planets, recent studies show that these 
models can also replicate the equatorial zonal flows of the four 
giant planets via the inclusion of an additional forcing mech- 
anism such as water vapor condensation (Lian and Showman, 
2010) or via enhanced radiative damping (Scott and Polvani, 
2008; Liu and Schneider, 2011). 

In the "deep-forcing" scenario, zonal winds are maintained 
by deep-seated convection and extend over the whole molec- 
ular envelope. These models rely on 3-D numerical simula- 
tions of rapidly-rotating spherical shells. In rotating convec- 
tion at moderate convective forcing (see Julien et al., 2012; 
King and Aurnou, 2012), convection occurs on long, axially- 
oriented columns reaching through the whole fluid layer. These 
well-organized columnar flows generate Reynolds stresses 
(i.e. statistical correlations between the convective flow com- 
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ponents) that drive strong zonal winds (e.g. Busse, 1983; 
Plautetal., 2008). The typical azimuthally prograde tilt of 
these convective columns always yields an eastward equatorial 
jet (e.g. Zhang, 1992; Christensen, 2001; Aurnou and Heimpel, 
2004). The direction and the number of jets are consistent 
with Jupiter's and Saturn's observation provided thin shells and 
low Ekman numbers are considered (e.g. Heimpel et al., 2005; 
Heimpel and Aurnou, 2007). 

A different zonal flow regime is however found when global- 
scale buoyancy starts to dominate the Coriolis force. The 
equatorial zonal flow tends to reverse (e.g. Gilman, 1977; 
Gilman and Foukal, 1979; Glatzmaier and Gilman, 1982). Re- 
sponsible for the latter is the turbulent mixing of angular mo- 
mentum, which may explain the strong retrograde equatorial 
zonal flow observed on the ice giants (Aurnou et al., 2007). 

To further inform the ongoing discussion on the driving 
mechanisms and the depth of the zonal jets in the giant plan- 
ets, here we investigate the deep-seated zonal flow perspective 
and thus compute 3-D global models of convection in spherical 
shells. 

Many of the previous parameter studies have employed the 
Boussinesq approximation where the density stratification is 
simply ignored (Christensen, 2002; Aurnou et al., 2007). This 
is rather dubious in the strongly stratified giant planets where 
the density contrasts are huge (e.g. Nettelmann et al., 2012a,b). 
More recent models therefore use the anelastic approxima- 
tion which allows to incorporate the effects of background 
stratification while filtering out the fast acoustic waves (e.g. 
Braginsky and Roberts, 1995; LantzandFan, 1999). In an 
extensive parameter study, Gastine and Wicht (2012) concen- 
trate on the influence of the density stratification on convec- 
tion and zonal flows in the rotation-dominated regime (see 
also Showman et al., 2011). While the density stratification 
affects the local scales and the amplitude of the convective 
flow, the mean zonal flows and the global quantities are fairly 
independent of the density contrast, similar to the results of 
Jones and Kuzanyan (2009). 

Many anelastic and fully compressible models of solar and 
stellar convection have observed a transition between the solar- 
like (i.e. eastward equatorial zonal flow) and the so-called 
"anti- solar" (i.e. westward equatorial zonal flow) differen- 
tial rotation profiles when buoyancy dominates the force bal- 
ance (Glatzmaier and Gilman, 1982; Bessolaz and Brun, 2011; 
Kapyla et al, 201 la; Matt et al., 201 1). Yet, to date, no system- 
atic parameter study has been made to investigate the influence 
of density stratification on the transition between the rotation- 
dominated and the buoyancy-dominated zonal flow regimes. 

This is precisely the focus of the present study, which ex- 
tends the previous Boussinesq study of Aurnou et al. (2007) to 
anelastic models. To this end, we conduct a systematic param- 
eter study from Boussinesq to strongly stratified models (i.e. 
Pbot/Ptop - 150) and solutions that span the range from weakly 
to strongly supercritical convection. 

In section 2, we present the anelastic formulation and the 
numerical methods. Section 3 shows the evolution of convec- 
tion when the driving is gradually increased from the rotation- 
dominated to the buoyancy-dominated regime. Section 4 fo- 



cuses on the zonal flows profiles that develop in the buoyancy- 
dominated regime. In section 5, we concentrate on the so- 
called transitional regime, a specific feature of strongly strat- 
ified anelastic models, before concluding in section 6. 



2. Hydrodynamical model and numerical methods 

2.1. Governing equations 

We consider hydrodynamical simulations of an anelastic 
ideal gas in a spherical shell rotating at a constant rotation rate 
Q about the z-axis. We use a dimensionless formulation of the 
governing Navier-Stokes equations where the shell thickness 
d - r - r t is employed as a reference lengthscale and Q" 1 as 
the time unit. Density and temperature are non-dimensionalised 
using their outer boundaries reference values p top and 7\ op . En- 
tropy is expressed in units of As, the imposed entropy contrast 
over the layer. Kinematic viscosity v, thermal diffusivity k and 
heat capacity c p are assumed to be homogeneous. 

Following the anelastic formulations of 
Gilman and Glatzmaier (1981); Braginsky and Roberts (1995) 
and Lantz and Fan (1999), the background reference state (de- 
noted with tildes in the following) is hydrostatic and adiabatic. 
It is defined by df/dr = -g/c p and a poly tropic gas p = f m , 
m being the poly tropic index. As we are interested in the 
dynamics of the molecular region of giant planets, we assume 
that the mass is concentrated in the inner part, so that g oc 1/r 2 
provides a good approximation (see also Jones and Kuzanyan, 
2009; Gastine and Wicht, 2012). Such a gravity profile then 
leads to the following background temperature profile 
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Here f and p are the background temperature and density, n = 
rijr is the aspect ratio of the spherical shell and N p corresponds 
to the number of density scale heights over the layer (see also 
Jones et al., 201 1, for the full derivation of the reference state). 

In the anelastic approximation, the dimensionless equations 
that govern convective motions are given by 



V • ipu) = 0, 



(3) 
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where w, p and s are velocity, pressure and entropy, respec- 
tively. S is the traceless rate-of- strain tensor with a constant 
kinematic viscosity, given by 
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6ij being the identity matrix. The dimensionless entropy equa- 
tion reads 



(ds \ E E 

— + u • Vs\ = - V • (pTVs) + — (1 - rj)c Q v , (6) 



where Q v is the viscous heating contribution given by 



2v = 2p 



(7) 



In addition to the aspect ratio n and the two parameters involved 
in the description of the reference state (N p and m), the system 
of equations (3-6) is governed by three dimensionless param- 
eters, namely, the Ekman number, the Prandtl number and the 
modified Rayleigh number: 
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where g top is the gravity at the outer boundary. Ra* can be re- 
lated to the standard Rayleigh number Ra = g top d 3 As/c p VK with 
(e.g. Christensen, 2002) 
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The definition of Ra* is based on the global entropy jump over 
the spherical shell and on the gravity value at the outer bound- 
ary. In anelastic models, it is however more appropriate to de- 
fine a local Rayleigh number that encompasses the radial de- 
pendence of the background reference state (e.g. Kaspi et al., 
2009; Gastine and Wicht, 2012): 



<R*(r) = 
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dr 
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where s c is the conduction state entropy, which is the solution 
of 



(pfVs c ) 



0. 
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As the entropy gradient is inversely proportional to pf , ft* can 
become very large at the surface in the most stratified cases. In 
these models, convection sets in first in the outermost region 
(Jones et al., 2009; Gastine and Wicht, 2012). To compare nu- 
merical models with different density contrasts, it is either pos- 
sible to consider a mass-weighted average of 7?* as suggested by 
Kaspi et al. (2009), or use its value at mid-depth (Unno et al., 
1960; Glatzmaier and Gilman, 1981). Table 1 shows that these 
two definitions lead to very similar results. In the following, we 
use the modified Rayleigh number at mid-depth 9^ = < ^*(r m i ( j) 
to compare our different numerical models. 

The first numerical models of rotating convection in spheri- 
cal shells by Gilman (1977) and Glatzmaier and Gilman (1982) 
have shown that the physical mechanism responsible of the 
zonal flow production is sensitive to the relative contribution of 
buoyancy and Coriolis force in the global- scale force balance. 
The ratio between these two forces can be related to Ra* (for 
the full derivation, see Aurnou et al., 2007) via 



Table 1: Local Rayleigh number at mid-depth and mass-weighted average of 
<R* for Ra* = 1 for different density stratifications. 
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Coriolis 

which is commonly referred to as the "convective Rossby 
number" in the solar and stellar convection communities (e.g. 
Elliott et al., 2000; Ballot et al., 2007) and in the fluid physics 
community (e.g. Zhong and Ahlers, 2010). In Boussinesq stud- 
ies, Ra* ~ 1 is typically found to be a good proxy to separate 
the rotation-dominated zonal flow regime (i.e. Ra* « 1) from 
the buoyancy-dominated flow regime (i.e. Ra* » 1) (Gilman, 
1977; Aurnou et al., 2007; Evonuk and Samuel, 2012). 

In contrast to the Ra* ~ 1 zonal flow transition in Boussi- 
nesq spherical shells, there is no consensus on the mecha- 
nisms that control the breakdown of the smaller- scale con- 
vection columns (e.g. Schmitz and Tilgner, 2009; King et al., 
2012; King and Aurnou, 2012; Julien et al., 2012). Carried out 
in Cartesian geometries in which zonal flows do not develop, 
these Boussinesq studies show that the columnar modes in ro- 
tating convection already tend to become unstable for values 
Ra* < 1 at low Ekman numbers. In our simulations, the convec- 
tion columns and the zonal flows seem to undergo simultane- 
ous behavioral transitions. This may occur due to the moderate 
Ekman values at which we carry out our suite of simulations, 
or due to a difference in the stability properties of convection 
columns in spherical shell convection in the presence of zonal 
flows. Thus, we will focus here only on the one fundamental 
behavioral transition that separates the rapidly-rotating regime, 
in which columns exist and the equatorial zonal flow tends to 
be prograde (e.g. Gastine and Wicht, 2012), and the buoyancy- 
dominated regime, in which the columns are unstable and the 
equatorial zonal flows tend to be retrograde (e.g. Aurnou et al., 
2007). 

2.2. Numerical methods and boundary conditions 

The numerical simulations of this parameter study have been 
carried out using the anelastic version of the code MagIC 
(Wicht, 2002; Gastine and Wicht, 2012), which has been val- 
idated in the Jones et al. (2011) anelastic dynamo benchmark 
study. To solve the system of equations (3-6), pu is decom- 
posed into a poloidal and a toroidal contribution 



pu = V x (V x We r ) + V x Ze r . 



(13) 



W, Z, s and p are then expanded in spherical harmonic 
functions up to degree £ max in colatitude and longitude 
(p and in Chebyshev polynomials up to degree N r in ra- 
dius. A detailed description of the complete numerical 
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Table 2: Critical Rayleigh numbers and corresponding azimuthal wave numbers 
for the different parameter regimes considered in this study. 
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3 x 1(T 4 
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1.291 x 10 5 
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3 x 10" 4 
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3.714 x 10 5 


47 


3 x 10" 4 


5 


6.373 x 10 5 


60 



method and the associated spectral transforms can be found 
in Gilman and Glatzmaier (1981). Typical numerical resolu- 
tions range from (N r = 65, £ max = 85) for Boussinesq models 
to (N r = 161, ^ max = 256) for the most demanding anelas- 
tic models with N p = 5. In the latter, a twofold or a fourfold 
symmetry in longitude has been used to ease numerical com- 
putation. As rapidly-rotating convection is dominated by high 
azimuthal wave numbers, this enforced symmetry is not consid- 
ered to be influential on the averaged properties of the flow (e.g. 
Christensen, 2002; Heimpel et al., 2005; Jones and Kuzanyan, 
2009). 

In all the numerical models presented in this study, we have 
assumed constant entropy and stress-free mechanical boundary 
conditions at both spherical shell boundaries, r t and r . 

2.3. A parameter study 

In this investigation, we consider two different Ekman num- 
bers, E = 10" 3 and E = 3 x 10" 4 , which are relatively large 
but allow us to carry out strongly supercritical cases. Follow- 
ing our previous hydrodynamical models, the Prandtl number 
is set to 1 and we use a polytropic index m = 2 for the reference 
state. For all the models of this study, we employ an aspect ratio 
rj = 0.6. Although this leads to deeper convective layers than 
the expected molecular to metallic hydrogen transitions in the 
gas giants (i.e. 0.85 Rj and 0.65 R s , see Liu et al., 2008), it has 
the advantage of avoiding the numerical difficulties associated 
with thinner shells. 

We have performed numerical simulations with various den- 
sity contrasts spanning the range from N p = 10" 2 (i.e. nearly 
Boussinesq) to N p = 5 (i.e. Pbot/Ptop - 150). The strongest 
stratification considered here is still below the expected density 
contrast of the gas giants interiors (i.e. N p ^ 8.5 from the 1 
bar level to the bottom of the molecular envelope at 1.5 Mbar, 
see Guillot, 1999; Nettelmann et al., 2012a). Since most of the 
density contrast is concentrated in the outermost region, a value 
of A/p = 5 already covers 99% of the radius of the molecular en- 
velope. For each density stratification and Ekman number, we 
vary the Rayleigh number from onset of convection to ~ 10. 
Critical Rayleigh have been obtained with the linear stability 
code by Jones et al. (2009) and are given in Tab. 2. 

Most of the runs have been initiated from a conductive ther- 
mal state (see Eq. 10) with a superimposed random entropy per- 



turbation. For the most demanding high ( FH m cases, a converged 
solution at different parameter values has been used as a start- 
ing condition. Altogether, more than 125 simulations have been 
performed, each running for at least 0.3 viscous diffusion time 
such that a nonlinearly saturated state is reached (see Tab. A.4). 

3. From rotation- to buoyancy-dominated regime 

3.1. Convective flows 

Figure 1 shows the evolution of the convective flow when the 
supercriticality is gradually increased for two different density 
stratifications. In the rotation-dominated regime (i.e. ^ « 1, 
upper panels of Fig. 1), the convective structures are aligned 
with the rotation axis following the Taylor-Proudman theo- 
rem. The azimuthal lengthscale of the convective columns is 
roughly three times smaller in the strongly stratified than in the 
Boussinesq model, following the critical wave numbers listed in 
Tab. 2. As shown in the linear stability analysis by Jones et al. 
(2009), this variation is due to the background density contrast 
that confines convection close to the outer boundary when N p 
increases. This confinement causes the decrease in the typical 
length scales in both the radial and the azimuthal directions (see 
also Gas tine and Wicht, 2012). 

Convection develops both inside and outside the tangent 
cylinder in cases with more supercritical Rayleigh numbers 
that still have <R* m < 1 (Fig. 1, middle row). The integrity of 
the convective columns is disturbed due to the gradual loss of 
geostrophy as the buoyancy forcing increases in strength (e.g. 
Soderlund et al, 2012). 

The alignment of the convective features along the rota- 
tion axis is completely lost when buoyancy dominates the 
global-scale force balance (^ > 1, lower panels of Fig. 1). 
While up- and downwellings have a very similar structure in 
the Boussinesq model, a strong asymmetry is visible in the 
strongly stratified case. Here the convection roughly forms a 
network of thin elongated downflows that enclose broader and 
weaker upflows. Due to the background density contrast, up- 
welling structures tend to expand and acquire a mushroom- 
like shape, while downwelling plumes are narrow and concen- 
trated. This network-like pattern of convection has been fre- 
quently observed in numerical models of the solar granulation 
(e.g. De Rosa et al., 2002; Miesch et al, 2008). 

3.2. Zonal winds 

Figure 2 shows the evolution of the zonal winds when < R^ n is 
increased in the nearly Boussinesq models (N p = 10" 2 , left pan- 
els) and in the strongly stratified models (N p = 5, right panels). 
The typical pattern of an eastward (i.e., prograde) outer and 
a westward inner geostrophic zonal flow already develops at 
mildly supercritical Rayleigh numbers (top row). These zonal 
flows are driven by Reynolds stresses (i.e. the statistical corre- 
lation of convective flow components), which rely on the pro- 
grade tilt of the convective columns induced by the boundary 
curvature (e.g. Busse, 1983; Zhang, 1992; Christensen, 2002). 
In the rotationally-dominated models, Reynolds stresses gener- 
ate a positive angular momentum flux away from the rotation 
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Figure 1: Radial velocity w r close to the outer boundary (r = 0.9 r Q ) for six different numerical simulations with E = 10 . Boussinesq models (i.e. with N p = 10 ) 
with increasing < R" m are displayed in the left panels, while anelastic models with N p = 5 are displayed in the right panels. Outward flows are rendered in red, inward 
flows in blue. 



axis that is balanced by viscous drag. Because of the strong 
confinement of convection close to the outer boundary in the 
N p = 5 case, the equatorial jet is narrower than in the Boussi- 
nesq model. In addition, its amplitude is significantly larger 
than the amplitude of the adjacent retrograde jet in the stratified 
model. 

With further increasing supercriticality, the discrepancies 
discussed above tend to vanish. Thus, for < R m < 0.5-1 
(third panels), the equatorial jets have indeed very similar am- 
plitudes and latitudinal extent while the retrograde jet is still 
somewhat weaker in the anelastic model. This agrees with the 
results or our previous parameter study in the < 1 regime 
(Gastine and Wicht, 2012). Note that the larger Ekman number 
cases considered here do not allow for the five jets solutions 
found in Gastine and Wicht (2012). 

When buoyancy starts to dominate the force balance (i.e. 
~ 1)' me zona l A° w direction largely reverses. The flow 
outside the tangent cylinder becomes mainly retrograde and 
the flow inside the tangent cylinder prograde. The amplitude 
of these zonal winds also increases roughly by a factor 5 at 



this transition reaching Ro ^ -0.5. The zonal winds are 
nearly z-independent for < R! m ^ 1 (fourth panels), even though 
the solutions are no longer in geostrophic balance. However, 
a further increase of the supercriticality produces some pro- 
nounced smaller-scale z-dependent features in the Boussinesq 
models. These small-scale structures are highly time-dependent 
and would therefore cancel out when time-averaged quantities 
are considered. The last panels of Fig 2 show that both the mean 
zonal flow and the small-scale ageostrophic motions depend on 
the density contrast in the » 1 regime. 

This evolution of the zonal flow structure is confirmed by the 
time-averaged toroidal kinetic energy spectra shown in Fig. 3 
for four Boussinesq models at E = 10" 3 . In the rotation- 
dominated regime (JC m < 1, upper left panel), the spectrum 
peaks at m = 0, illustrating the strong contribution of ax- 
isymmetric zonal flows in the toroidal kinetic energy bud- 
get. For spherical harmonic orders m > 20, the spectra fol- 
lows a clear power-law behaviour and the steep decrease (close 
to m~ 5 ) is very similar to previous quasi-geostrophic studies 
(Schaeffer and Cardin, 2006) and matches the k~ 5 scaling de- 
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Figure 2: Zonally averaged azimuthal velocity in the meridian plane for simulations with E = 10~ 3 and increasing K^. Upper panels correspond to Boussinesq 
models (i.e. N p = 10~ 2 ) and lower panels to N p = 5. Colorscales are centered around zero: prograde jets are rendered in red, retrograde jets in blue. In some cases, 
the prograde contours have been truncated in amplitude to emphasise the structure of the retrograde flows. Extrema of the zonal flow velocity are indicated in the 
center of each panel (velocities are expressed in Rossby number units, i.e. u/£lr ). 
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rived by Rhines (1975) in the /3-plane turbulence framework. 
Axisymmetric zonal flows still dominate the toroidal energy in 
the K m = 2.92 and K m = 10.23 cases, but the slope at smaller 
lengthscales flattens gradually. At = 29.24 (lower right 
panel), the spectrum shows a mild maximum at m = 2 before 
clearly following an inertial range scaling of m" 5/3 , the theo- 
retical behaviour expected for homogeneous and isotropic 3-D 
turbulence (e.g. Lesieur, 1997). In this case, the axisymmet- 
ric zonal flow contribution is not dominant anymore as already 
suggested by the last panel of Fig. 2a. 

3.3. Convection regimes 

To further investigate the regime change observed in Figs. 1- 
2, we consider time-averaged quantities for the whole set of 
numerical models computed in this parameter study. We fo- 
cus in the following on the zonal flow amplitude characterised 
by the dimensionless Rossby number Ro = /Qr and on the 
contribution of axisymmetric toroidal energy in the total kinetic 
energy budget. Following Christensen (2002), this is quantified 
by the ratio of the total to the non-zonal kinetic energy E^ n /E nz , 
where £kin has been obtained by 



Figure 3: Time-averaged toroidal kinetic energy spectra for four different nu- 
merical models with N p = 10 -2 and E = 10 -3 . 
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Figure 4: (a) Amplitude of the surface zonal wind at the equator (in units of 
Ro e = U(f,/Qr ) plotted against K^. (b) Ratio of total over non-zonal kinetic 
energy plotted against 9^. In both panels, closed symbols correspond to simu- 
lations with E = 3 x 10~ 4 and open symbols to simulations with E = 10 -3 . The 
transition at ^ = 1 is emphasised by a dotted vertical line. 



where V is the volume of the spherical shell and r is the time- 
averaging interval. Each numerical simulation has been aver- 
aged long enough to suppress the short term variations. 

Figure 4 shows how the surface equatorial zonal flow ampli- 
tude Ro e and the ratio E^ n /E nz evolve when < R? m is increased 
for numerical models with different Ekman numbers and den- 
sity contrasts. In the rotation-dominated regime < 1), the 
equatorial jets is always prograde and its amplitude gradually 
increases with K^. A maximum amplitude of Ro e ^0.1-0.15 
is reached around ^ 0.5-0.8. In this regime, the influ- 
ence of the density contrast is not obvious as, on the one hand, 
the N p = 3 cases tend to produce stronger jets than the Boussi- 
nesq models, which are, on the other hand, rather similar to the 
strongly stratified N p = 5 cases. 

Figure 4b shows that the ratio E^ n /E nz first increases for 
weak to moderate supercriticalities before decaying once con- 
vection is strongly driven. This decay is attributed to the grad- 
ual decorrelation of the convective flow components in the 
strongly-driven regime, which reduces the efficiency of the en- 
ergy transfer between small-scale convection and large-scale 
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Figure 5: Regime diagram of the evolution of the mean surface zonal flow at the 
equator as a function of 9^. Regime I corresponds to the rotation-dominated 
regime, in which ^ < 1 . Regime II corresponds to the buoyancy-dominated 
regime in which the mixing of the angular momentum per unit of mass is ef- 
ficient, while Regime III suggests a possible suppression of the zonal flow for 
ni » 1. 



mean zonal flows via Reynolds stresses (Christensen, 2002; 
Gastine and Wicht, 2012). The maximum of this ratio depends 
on both the density stratification and the Ekman number. In 
Boussinesq models at E = 3 x 10" 4 , this maximum is around 10 
(which means to 90% of the total kinetic energy is contained in 
the zonal flows) while it reaches only 5 (i.e. 80%) in the N p = 5 
models. In addition, decreasing the Ekman number tends to 
produce a stronger zonal flow contribution (see Christensen, 
2002). 

An abrupt transition to retrograde equatorial zonal winds 
takes place close to ^ 1 , independently of the density strat- 
ification and the Ekman number considered. The retrograde 
equatorial jet amplitude is roughly multiplied by a factor 4 at 
this transition to reach Ro e ^ -0.4. This goes along with 
a larger E^ n /E nz value that reaches approximately 7 for the 
E = 10" 3 cases (i.e. 85%) and 15 for the E = 3 x 10" 4 cases 
(i.e. 93%). 

As is increased further, the density stratification starts 
to play a more important role. For stronger stratifications, 
the equatorial zonal flow amplitude further increases while it 
roughly levels off and finally decreases in the Boussinesq mod- 
els. This is also reflected in the evolution of Ey± n /E nz where 
the strongly stratified cases decay slower at large than the 
weakly stratified models. Since we expect that the convective 
fluctuations will dominate the mean zonal flows at ^ » 1 , the 
ratio £kin - E nz should ultimately tend to unity. This is already 
observed in Boussinesq cases in Fig. 4b where a value of 1 .04 is 
reached at = 29.24. For the stronger stratified cases, much 
larger < R^ n values are required than we could afford to simulate 
numerically. 

Figure 5 shows a tentative regime diagram based on our sim- 
ulation results. Regime I corresponds to rotation-dominated 
cases i^RL < 1) where convection shows a columnar structure 
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and maintains a prograde equatorial zonal flow. Regime II is 
characterised by three-dimensional convection and a strong ret- 
rograde equatorial jet. The transition between regimes I and II 
takes place around < K m ~ 1 independently of the density stratifi- 
cation and the Ekman number. The strong decrease of the zonal 
flow amplitude observed in the Boussinesq models for ( K m > 5 
points toward a possible third regime in which the zonal flows 
may be gradually suppressed by the effects of turbulent diffu- 
sion. 

The reason for the density-dependent jets amplitude observed 
in the buoyancy-dominated regime and the possibility of a third 
regime at large Ra* are further discussed in the next section. 

4. Angular momentum mixing 

4.1. Influence of the density stratification 

When buoyancy dominates Coriolis force, strong convection 
can locally mix the angular momentum. Flows in which angular 
momentum is indeed homogenised have been found in several 
studies of atmospheric dynamics or numerical models of rotat- 
ing convection with a weak rotational influence (e.g. Gilman, 
1977; Hathaway, 1982; DeRosaetal., 2002; Aurnouetal., 
2007; Bessolaz and Brun, 2011). 

In an inertial reference frame, the angular momentum of a 
fluid parcel is given by 

M = Mzf + Mpr = pu<ps d^V + pas 2 d^V, (15) 

where S^V is the volume of the fluid element, s = rsinO is 
the cylindrical radius, Mzf is the angular momentum due to 
zonal flows and AIpr is the angular momentum that comes from 
the planetary rotation. The overbars here indicate an azimuthal 
average. Because of the anelastic continuity equation (3), the 
mass of a fluid element is conserved during its displacement 
such that pd'V = const.. If, in addition, a fluid parcel con- 
serves its angular momentum M, then, as hypothesized by 
Gilman and Foukal (1979), the angular momentum per unit of 
mass £ is a conserved quantity 



Table 3: £ from Eq. (18) and frc for different density stratifications. 
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Fi gure 6: Evolution of £ for different density stratifications as a function of s m i x > 
the minimum cylindrical radius above which mixing of angular momentum per 
unit of mass occurs. The vertical line corresponds to the location of the tangent 
cylinder. 



and decreases by 25% when the density contrast is increased 
from nearly Boussinesq (N p = 10" 2 ) to N p = 5 cases (see the 
middle column of Tab. 3). This decrease of £ explains the dif- 
ferences in equatorial jet amplitudes in Fig. 4 via 



r s 
Ro = C— . 

s r a 



(19) 



Note that Aurnou et al. (2007) derived the Boussinesq version 
of Eq. (18), in which the background density is homogeneous. 



L - Xzf + Xpr = u<ps + Qs = const. (16) 
Nondimensionalising this equation by £lr 2 leads to 

X* =Ro- + 4- (17) 

r o n 

A spatially homogeneous X* value would then equal the mass 
integral of the initial angular momentum distribution (i.e. a 
rigid body rotation) given by 



f = <r>, 



= li £ ™ dm 

m Jo Jo Jr, 



\p(r)s 2 



(18) 



r sin dr d6 dcf), 



where m is the total mass of the spherical shell given by m = 
J y p(r)dV. £ depends on the background density stratification 



4.2. Partial mixing of£* 

The previous derivation provides an idealised description of 
the simulation results at » 1. For instance, the equato- 
rial zonal flows, shown in Fig. 4, never reach the maximum 
theoretical amplitudes of Ro e - -0.528 for N p = 10" 2 and 
Ro e = -0.654 for N p = 5 predicted by Eq. (19). From this 
we infer that the mixing of £* takes place in only part of the 
spherical shell. 

We can then make the following ad-hoc assumption that X* 
is homogenised for cylindrical radii s > s m { x only. This leads to 



=— f f r 

fornix) J z(s) Jo J Sai 



\p(r)s 2 



sdsd(pdz, (20) 



where z(s) = ± ^r 2 - s 2 outside the tangent cylinder and 
± ^ y/r% - s 2 - ^Jr 2 - s 2 j inside. Figure 6 shows the evolution 




Figure 7: Left: time-average angular momentum per unit of mass plotted against the cylindrical radius s for various numerical simulations with E = 10~ 3 . Right: 
corresponding azimuthally-averaged velocity profiles as a function of latitude at the outer boundary. The vertical lines in the left panels correspond to the tangent 
cylinder. For comparison, the grey-shaded area are delimited by the theoretical values associated with £* mixing over the whole spherical shell (lower bound) and 
outside the tangent cylinder only (upper bound, see Tab. 3 for the corresponding values). 



of £ Smix for different density stratifications as a function of s m { x . 
When the shell is only partially mixed the influence of the 
density stratification is gradually reduced and becomes negli- 
gible when £,* is homogenised in only a small fraction of the 
domain. Because the tangent cylinder is often considered to 
present a dynamical barrier for rotation-dominated convection, 
Aurnou et al. (2007) introduced the respective values for mix- 
ing outside the tangent cylinder only. Note, however, that this 



is a rather arbitrary choice since there is no expected sharp vor- 
ticity step across the tangent cylinder in the regime II discussed 
here. The second column in Tab. 3 corresponds to these cases 
(also see the dotted lines in Fig. 6). In the Boussinesq limit we 
recover £ T c = l-f(l - ^ 2 ) already derived by Aurnou et al. 
(2007). 

Figure 7 displays time-averaged radial profiles of £* and lat- 
itudinal profiles of surface zonal flows for numerical simula- 
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C pu m p u m • VC = - V • J* tot = - V • F Re + V • J* vi[ 
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Figure 8: Azimuthal force balance for two numerical models with ^ = 0.24 (upper panels) and K* m = 2.37 (lower panels). Both cases have N p = 5 and E = 10 -3 . 
(a) Time-averaged angular momentum per unit of mass £,*. (b) Time-averaged stream-function of the meridional circulation pu m . (c) Advection of the angular 
momentum £* by the mean meridional circulation, i.e. left-hand side of Eq. (21). (d) The net axial torque of the right-hand side of Eq. (21). (e) The contribution of 
Reynolds stresses to this torque, (f) The contribution of viscous stresses to this axial torque. For each panel, positive (negative) values are rendered in red (blue). 
Panels (d-f) share the same contour levels. 



tions in the buoyancy-dominated regime (Regime II). For each 
panel, the grey-shaded area show the spans between the theo- 
retical X* -mixing over the entire shell and outside the tangent 
cylinder only. For comparison, a solid body rotation profile is 
demarcated by the dashed black lines in the left panels of Fig. 7. 

For < F? m *t 1 (dark green lines), the X* profiles are roughly 
constant outside the tangent cylinder in the Boussinesq models. 
In contrast, X* -mixing occurs only in a relatively thin region 
(s > 0.8 r ) in the strongly stratified cases. The Ro e ~ -0.4 
results are consistent with our predictions based on on Fig. 6 
(i.e. Ro e = £ Smix - 1) for the range of mixing depths that are 
observed. 

In the strongly stratified cases (N p > 3), the X* -mixing 
moves deeper when ( Ft m is increased. The entire region outside 
the tangent cylinder is homogenized (light green line for N p = 3 
and orange line for N p = 5) and the surface zonal flows match 
the theoretical profile derived before, at least at mid latitudes 
(i.e. 6 < ±50°) once K m ^ 2. For K m > 10, £* is constant 
for s > 0.5 r and the zonal flow profiles lie roughly midway 
between the two theoretical curves (fully-mixed and s m { x = stc 
limits). They reach stronger amplitude in the N p = 5 cases due 



to the larger available X* reservoir (Fig. 6). 

The profiles also show that the X* -mixing stops in the highest 
< R^ n Boussinesq models (dark purple line) and gradually tends 
toward the solid body rotation , i.e. £* = s 2 lr 2 . As described 
before, this weakening of the zonal flow at the highest accessi- 
ble values of <R* m appears to define a third dynamical regime. 

4.3. Azimuthal force balance and meridional flow structure 

In this section, we explore how the zonal flows are main- 
tained in the X* -mixing regime (Regime II in Fig. 5) and the 
possible decay of these flows at larger (Regime III). To do 
so, we analyse the axisymmetric, azimuthal component of the 
Navier-Stokes equations (Eq. 4): 



P— +pu m • VX 
ot 



= -V • T t ou 

= -v • tr Re - r V i SC ] , 

psu'u' - Eps 2 Vl — ] 



(21) 



= -V- 
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where u m is the mean meridional circulation, L = su^ + s 2 
is the angular momentum per unit of mass, Ttot is the angular 
momentum flux associated with Reynolds and viscous stresses, 
Trq and Tvisc, respectively. Primed quantities correspond to 
fluctuations about the axisymmetric mean. Equation (21) re- 
quires that the advection of axisymmetric angular momentum 
by the meridional flow is balanced, on time-average, by the net 
axial torques due to Reynolds stresses and viscosity. 

Figure 8 shows the angular momentum per unit of mass and 
the mean meridional flow, along with the different contributions 
of this force balance for two strongly stratified numerical mod- 
els with R* m = 0.24 and <R* m = 2.37. In the rotation-dominated 
case (*R* m = 2.37, upper panels), the equatorial zonal wind is 
prograde and L gradually increases with cylindrical radius s. 
In the buoyancy-dominated case = 2.37, lower panels), L 
is roughly constant outside the tangent cylinder (see also the 
orange curve in the last row of Fig. 7) and is associated with a 
strong retrograde equatorial jet which reaches Ro e ^ -0.5. 

The meridional circulation patterns differ between the two 
numerical models shown in Figure 8. While a relatively small- 
scale multicellular meridional circulation structure is observed 
in the first case, the meridional flow in the second model is 
dominated by a pair of large-scale cells. In the latter, a second 
pair of weaker counter-cells is discernable close to the inner 
boundary. Such a transition between small-scale multicellular 
and large-scale single-celled meridional circulation has been al- 
ready observed for numerical simulations around < Ft m ~ 1 (e.g. 
Elliott et al, 2000; Matt et al, 201 1; Bessolaz and Brun, 201 1). 

This change of the meridional circulation pattern observed 
around ^ 1 is reflected in the spatial variations of the force 
balance expressed by Eq. (21). The third and fourth panels of 
Fig. 8 check this balance and confirm that the advection of L by 
the mean meridional flow is indeed balanced, on time average, 
by the sum of the torques due to viscous and Reynolds stresses. 
Some small-scale time-dependent features remain visible in the 
= 2.37 case but will vanish if a longer averaging time is 
considered. In the = 0.24 case, viscous stresses (panel f) 
are significant and largely, but not perfectly, compensate the 
Reynolds stresses (panel e). This results in several sign changes 
in the axial torque (panel d) that are mirrored by the com- 
plex meridional flow structure (see also Augustson et al., 2012). 
This force balance is somewhat different in the - 2.37 
model in which the simple large-scale Reynolds stresses domi- 
nate the force balance and drive the dominant pair of meridional 
circulation cells. The viscous contribution becomes significant 
only near the inner boundary and results in the second weaker 
pair of cells visible in Fig. 8b. 

4.4. Towards a possible third regime for < R^ n » 1 

Figure 8 shows that the Reynolds stresses become the dom- 
inant contribution to the net axial torque when < R? m > 1 . The 
force balance (21) can thus be approximated in the < K m » 1 
regime by 




Figure 9: Correlation C r ^> for different numerical models in regime II with N p = 
10~ 2 , E = 10 -3 (upper panels) and N p = 5, E = 10 -3 (lower panels). The 
average correlation C r( p throughout the shell is indicated in the middle of each 
panel. 



Convection Meridional circulation V~7~\ Zonal flows 




ft* = 1.14 ft* = 1.76 ft* = 3.52 ft* = 6.16 ft* = 17.59 




VX - -V • r Re - -V • [ps u'u'^ , (22) 



ft*=0.71 ft* =1.19 ft* =2.37 ft* =7.12 ft* = 11.87 



Figure 10: Time-averaged kinetic energy budget for the numerical simulations 
displayed in Fig. 7. Dark- grey area correspond to the non-zonal kinetic en- 
ergy, black area to the kinetic energy contained in the axisymmetric poloidal 
flows (i.e. meridional circulation), and hatched light grey area to the energy in 
axisymmetric toroidal flow. 
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since viscosity does not play a significant role here. This 
balance shows that Reynolds stresses rely on the correla- 
tions between the meridional and the longitudinal components 
of the convective flow u' r u'^ and u' B u'^ (e.g. Ruediger, 1989; 
Kapyla et al., 2011b). In the following, we focus on the evo- 
lution of u' r u' only and quantify the correlation by 



U r U '(j) 

C ^ = /__ x l/2' (23) 
(< 2 <) 

where the overbars correspond to an azimuthal average. Fig- 
ure 9 shows the evolution of C r( p when < R$ n is increased in 
Boussinesq (upper panels) and in strongly stratified models 
(lower panels). For < R m 1, the correlations are significant and 
strong negative Reynolds stresses are maintained. However, C r( p 
gradually decreases in the more supercritical cases. An increase 
of the Rayleigh number indeed goes along with stronger tur- 
bulent velocities and smaller typical flow lengthscales, leading 
to a gradual decrease of the turnover timescale of convection. 
Small-scale eddies are not influenced by rotation anymore since 
their lifetime becomes significantly smaller than the rotation 
period (see also Gastine and Wicht, 2012). This loss of coher- 
ence results in a gradual decrease of the Reynolds stresses cor- 
relations needed to maintain the mean flows (Brummell et al., 
1998; Mieschetal., 2000; Kapyla etal., 2011b). A reduction 
of the mean flows amplitude is therefore anticipated when the 
contribution of the turbulent motions increases in the total en- 
ergy budget. 

Figure 10 shows the energy distribution for numerical mod- 
els in the buoyancy-dominated regime (see Fig. 7 for the corre- 
sponding zonal flows). Zonal flows dominate the energy budget 
for ( K m ^ 1 (80-90%). In Boussinesq models, the contribution 
of the turbulent flows then increases rapidly to overwhelm the 
energy distribution for > 10. This increase is more gradual 
for the N p = 3 cases and it even seems to level around 20% 
in the N p = 5 models. This confirms the density-dependent 
evolution of Ey^ n /E nz in the » 1 regime visible in Fig. 4b. 

This gradual decorrelation defines a transition towards the 
possible regime III displayed in Fig. 5. When the degree of tur- 
bulence increases and the small-scale motions dominate the en- 
ergy budget, a decrease of zonal flows amplitudes is observed in 
Boussinesq and N p = I models due to the gradual loss of corre- 
lation of the convective flow. However, this scenario still needs 
to be confirmed in strongly stratified models: the presently ac- 
cessible range of < R m allows for only weakly turbulent convec- 
tive motions in our N p > 3 simulations (see Tab. A.4). 



5. A transitional regime in anelastic models 

As mentioned in section 2, <R* is a radially-dependent quan- 
tity that can vary across the fluid shell by several orders of mag- 
nitude in the strongly stratified models. This can lead to differ- 
ent dynamical regimes close to the outer boundary and in the 
fluid shell's deeper interior. 



io- 1 





^^TC > 1 


TV <\^^^ 


dx a) 



0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 




0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 



r/r 

Figure 11: (a) Sketch of the radial profile of K* in strongly stratified models. 
r mix corresponds to the radius at which <R* = 1. It marks the tentative limit 
between the buoyancy-dominated region (<R* > 1) and the rotation-dominated 
region (<R* < 1), both emphasised by a grey-shaded area, (b) Radial profile of 
<R* for various numerical simulations with N p = 5 and E = 3 x 10 -4 . 




Figure 12: Sketch illustrating the transitional regime in strongly stratified mod- 
els. 
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Figure 13: Vorticity along the axis of rotation co z = (V x u) z displayed in the equatorial plane for three numerical models with N p = 5 and E = 3 x 10 4 . Red (blue) 
correspond to positive (negative) values. The dashed black circles correspond to values defined in Fig. 11. 



5.7. Convective flows in the transitional regime 

Figure 1 1 shows the evolution of *R* (r) for different numeri- 
cal models with N p = 5. The spherical shell can be separated in 
two distinct layers: > 1 close to the outer boundary where 
buoyancy effects become large and <R* < 1 in the deep inte- 
rior where rotation can still dominate the force balance. As 
( R" m increases, the buoyancy-dominated region grows inward. 
The radius r m i x is defined by the radius at which <R* = 1 ap- 
proximately separating the buoyancy-dominated region from 
the rotation-dominated inner region. Figure 12 sketchs this 
transitional regime in which columnar convection in the deep 
interior exists contemporaneously with three-dimensional con- 
vective structures close to the outer boundary. 

Figure 13, showing equatorial cuts of oj z for numerical mod- 
els with N p = 5 and three different Rayleigh numbers, confirms 
that r m i x indeed coincides with a dynamical regime bound- 
ary. In the inner part, the convective flow is dominated by 
convective columns tilted in the prograde direction (positive 
vorticity dominates) that maintain Reynolds stresses. Beyond 
r m i x , the convective flow seems to be more radially-oriented 
with no preferred sign of vorticity anymore. Some strongly- 
stratified simulations discussed in our previous study already 
show the very beginning of this transitional regime (see Fig. 14 
in Gastine and Wicht, 2012). 

5.2. Zonal flows in the transitional regime 

The transition to regime II happens when r m i x reaches ap- 
proximately mid-depth, where an abrupt transition to a retro- 
grade equatorial jet takes place. The parameter space, where 
the transitional regime with two distinct types of convection co- 
exists, is therefore quite narrow. Strong stratification is required 
and has to be neither too small nor too large, typically in the 
K m ^ 0.1-0.5 range. 

Figure 14 shows the surface zonal flows for numerical mod- 
els that lie in this specific range of parameters. Typical for nu- 
merical simulations at relatively large Ekman numbers (here 
E = 3 x 10" 4 ) that are still in the < 1 regime, a large 
prograde equatorial jet is flanked by two weaker retrograde 
jets at higher latitudes (±60°). Due to the aspect ratio of the 
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Figure 14: Time-averaged surface zonal flows as a function of latitude for var- 
ious numerical simulations with N p = 5 and E = 3 x 10~ 4 . The corresponding 
K*{r) profiles are given in Fig. 11. 



spherical shell considered in this study (tj = 0.6), the pro- 
grade equatorial wind extends up to mid latitudes of ±50°. For 
< R0 m - 0.128 (green line), the zonal wind maximum is reached 
at the equator, similarly to the typical Boussinesq models (e.g. 
Heimpel et al., 2005). However, when increasing K^, the cen- 
ter of the main equatorial jet decreases until eventually a dim- 
ple appears flanked by two maxima. The width of the dim- 
ple further grows with < R m while the central flow amplitude de- 
creases reaching approximately ±25° latitude and 20% of the 
maximum zonal wind amplitude just before the transition to 
regime II (see the inset in Fig. 14). Such a pronounced dimple 
has also been observed in the anelastic models by Kaspi et al. 
(2009) who also employed strong density contrasts (see Fig. 13 
in their study). Following the sketch displayed in Fig. 12, this 
dimple can be directly related to the transition between rotation- 
dominated and buoyancy-dominated regions. As already shown 
in Fig. 13, r m i x separates these two regions and thus allows us 
to roughly estimate the typical latitudinal extent of the dimple 
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Figure 15: Observed surface zonal winds on Jupiter. Velocities are given 
in Rossby number units. The data are adapted from Porco et al. (2003) and 
Vasavada and Showman (2005). 

via a simple geometrical relation (see Fig. 12) 

COS ^dimple • (24) 

r o 

This simple expression thus relates the latitudinal extent of the 
observed surface zonal flows to the transition radius between 
the two distinct internal dynamical regimes. 

5.3. A dimple in the equatorial zonal band of Jupiter 

A similar dimple exists on Jupiter (Fig. 15) and extends 
roughly between ±7° latitude. The amplitude of the equato- 
rial wind decreases by roughly 30% at the equator. Is Jupiter's 
dimple the surface expression of an internal regime transition? 
From Eq. (24) with Simple ~ 7°, we can speculate that the mix- 
ing radius r m i x would lie around 500 km below the 1 bar level. 

To evaluate the plausibility of this scenario, we try to es- 
tablish the value 7?* urf of Jupiter using the scaling laws by 
Christensen (2002) and Christensen and Aubert (2006) l , which 
relates the modified Nusselt number to the flux-based Rayleigh 
number Raj = agq/pc p Q?d 2 via Nu* = 0.076 Raj 053 . 7?* urf 
can then directly be obtained with 7?* urf = Ra*/Nu*. Tak- 
ing an internal heat flux of 5.5 W.m -2 , a = 6 x 10~ 3 K _1 , 
Q = 1.75 x lO-V 1 , c p = 1.2 x lO^.kg^K" 1 , a gravity 
g = 25m.s~ 2 , a density p = 0.2kg.m~ 3 and a lengthscale 
d = 0.1r jup = 7 x 10 6 m (values taken from Guillot, 1999; 
French et al., 2012), we derive a surface value 7?* urf ~ 10" 2 . 

As suggested by Miesch and Hindman (2011), another way 
to evaluate the impact of rotation on convection from observ- 
able quantities is to measure the convective Rossby number 
at the surface Ro c = (Ot c ) _1 , where r c is the convective 
turnover time (see also Gastine and Wicht, 2012, for some com- 
ments on this quantity as a viable proxy to quantify the tran- 
sition between the two regimes). Using r c = H p /u com with 



! Note that using Nu ~ Ra 1/3 as suggested by King et al. (2012) leads to 
very similar values of Ra*. 



H p ~ 5 x 10 4 m from the models by Nettelmann et al. (2012a) 
and u com ~ 5 m.s -1 from observations by Salyk et al. (2006), 
we obtain r c ~ 10 4 s, which leads to Ro c ~ 0.5 (i.e. 7?* urf ~ 
(Ro c ) 2 ~ 0.25) at Jupiter's surface, an order of magnitude larger 
than the scaling law prediction. 

Concerning the uncertainties in our simulations as well as in 
our knowledge of Jupiter's dynamics, this latter value is already 
quite close to the expected value of ^* urf > 1 required to form a 
dimple in the tansitional regime. 

6. Conclusion 

We have investigated the transition between the rotation- 
dominated and the buoyancy-dominated regimes in rotating 
spherical shells with different density contrasts, extending the 
previous Boussinesq study by Aurnou et al. (2007). Follow- 
ing Gilman and Glatzmaier (1981) and Jones and Kuzany an 
(2009), we have employed the anelastic approximation to filter 
out fast acoustic waves and the related short time steps. Explor- 
ing moderately small Ekman numbers (E = 10" 3 - 3 x 10" 4 ) 
allowed us to raise the Rayleigh number into the buoyancy- 
dominated regime (characterised here by < K m > 1) and to study 
zonal flows in a broad parameter range. We highlight our main 
findings: 

• When gradually increasing K^ 9 the convective flows 
evolve from geostrophic columnar convection when rota- 
tion dominates the force balance to three-dimensional tur- 
bulent motions when buoyancy becomes significant. In the 
stratified cases, the latter is characterised by a pronounced 
asymmetry between broad upwellings and narrow down- 
wellings. This evolution of the convective features is ac- 
companied by a sharp transition in the zonal flow regime. 
The equatorial zonal jet reverses its direction at < R? m ^ 1, 
independently of the background density contrast and the 
Ekman number. 

• In the rotation-dominated regime (i.e. «: 1), a com- 
bination of geostrophic columns and boundary curvature 
lead to Reynolds stresses (i.e. the correlation between the 
cylindrically radial and the azimuthal components of the 
convective flow) that maintain a prograde equatorial zonal 
flow. The zonal flow amplitude is largely independent 
of the density stratification (see also Gastine and Wicht, 
2012). 

• In the buoyancy-dominated regime (i.e. { K m » 1), con- 
vection homogenises the angular momentum per unit of 
mass which leads to a strong retrograde equatorial zonal 
flow flanked by prograde winds at higher latitudes. The 
maximum zonal flow amplitude now increases with den- 
sity stratification. As already mentioned by Aurnou et al. 
(2007), these zonal flow patterns are reminiscent to those 
observed on Uranus and Neptune, though it remains un- 
certain whether convective driving in the ice giants indeed 
reaches K m > \. 



14 



• Our simulations suggest the possible existence of a third 
regime where the zonal flow amplitude decreases and three 
dimensional turbulent convection strongly dominates. The 
timescale of the small-scale convective motions is much 
shorter than the rotation period which therefore cease to 
play a role in organising large-scale flow. The transition to 
the third regime seems to depend on the density stratifica- 
tion and has not been reached for N p > 3. 

• For strongly stratified models in the < R m ^ 0.1 - 0.5 range, 
both the dynamical regimes I and II can be present in 
the spherical shell. Close to the outer boundary, buoy- 
ancy dominates and leads to more turbulent flows within a 
thin outer layer. In the deep interior, rotation still dom- 
inates, and columnar convection drive the typical zonal 
flow structure. The turbulent outer layer reduces the zonal 
flow amplitude in the center of the equatorial jet. This 
leads to a dimple similar to the one observed on Jupiter. 
Its width suggests that a transition between the two dy- 
namical regimes may occur at a depth of 500 km below 
the 1 bar level in Jupiter. Estimate based on the turnover 
timescale of convection suggests 7?* urf ~ 0.25, just on the 
low side of the required value. 

• This parameter study on the different zonal flows regimes 
in rotating anelastic spherical shells has also some possi- 
ble stellar applications as it helps to clarify the transition 
between the solar-like differential rotation (i.e. prograde) 
and the so-called anti-solar differential rotation (i.e. ret- 
rograde) observed in many previous stellar models (e.g. 
Bessolaz and Brun, 201 1; Kapyla et al., 201 la; Matt et al., 
2011). 

Theoretical and laboratory studies in Cartesian geome- 
tries find that the transition between rotation-dominated and 
buoyancy-dominated local-scale convection occurs at Ra* <?c 1 
as the Ekman number is made small (e.g. Julien et al., 2012; 
King et al., 2012). However, it is not clear how these findings 
relate to the system- scale zonal flow transitions investigated 
herein. Thus, further work is required to bridge the gap between 
the low Ekman number, Boussinesq, Cartesian studies and the 
moderate Ekman number, anelastic, spherical shell simulations 
we consider here. 
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Table A.4: Results table. 



Ekman 




Kn 


Ra/Ra c 


R^zon 


R^mer 


Re' 


Ro e 


Nu 


10~ 3 


O.Ol 


2.92 x 10~ 2 


1.4 


3.2 


0.0 


5.3 


1.85 x 10~ 3 


1.06 


10- 3 


0.01 


4.39 x 10~ 2 


2.1 


7.0 


0.1 


7.7 


3.64 x 10~ 3 


1.11 


1(T 3 


0.01 


7.31 x 10~ 2 


3.5 


14.7 


0.3 


12.6 


6.84 x 10~ 3 


1.21 


10~ 3 


0.01 


1.17 x 10 _1 


5.7 


29.3 


0.8 


20.8 


1.62 x 10~ 2 


1.44 


10" 3 


0.01 


1.46 x 10 _1 


7.1 


43.4 


1.1 


26.4 


2.55 x 10~ 2 


1.60 


10~ 3 


0.01 


2.19 x 10 _1 


10.6 


77.6 


3.0 


41.6 


4.59 x 10~ 2 


2.18 


10- 3 


0.01 


2.92 x 10 _1 


14.2 


107.5 


5.3 


57.6 


6.32 x 10~ 2 


2.88 


10~ 3 


0.01 


4.39 x 10 _1 


21.2 


145.4 


10.7 


90.0 


8.32 x 10~ 2 


4.22 


10- 3 


0.01 


5.85 x 10 _1 


28.3 


150.9 


11.6 


121.1 


1.05 x 10 _1 


5.36 


10~ 3 


0.01 


7.31 x 10 _1 


35.4 


154.1 


14.9 


149.6 


1.08 x 10 _1 


6.46 


10- 3 


0.01 


8.77 x 10 _1 


42.5 


140.3 


18.1 


178.7 


1.05 x 10 _1 


7.48 


10" 3 


0.01 


1.02 


49.6 


127.8 


21.0 


205.5 


1.01 x 10" 1 


8.37 


10" 3 


0.01 


1.17 


56.6 


511.1 


51.5 


201.5 


-3.73 x 10" 1 


8.32 


10- 3 


0.01 


1.32 


63.7 


538.2 


57.1 


215.4 


-3.82 x 10- 1 


9.16 


10" 3 


0.01 


1.46 


70.8 


561.3 


58.6 


231.2 


-3.85 x 10" 1 


9.90 


10- 3 


0.01 


2.92 


141.6 


614.2 


94.6 


376.6 


-4.08 x 10" 1 


15.10 


10" 3 


0.01 


4.39 


212.4 


548.8 


114.5 


483.6 


-4.06 x 10" 1 


18.70 


10- 3 


0.01 


5.85 


283.2 


471.6 


142.8 


583.5 


-3.88 x 10" 1 


21.70 


10" 3 


0.01 


7.31 


354.0 


460.6 


153.3 


650.5 


-3.57 x 10" 1 


23.00 


10- 3 


0.01 


1.02 x 10 1 


495.7 


432.1 


191.6 


802.9 


-3.03 x 10" 1 


27.00 


10- 3 


0.01 


1.46 x 10 1 


708.1 


351.7 


220.9 


1010.5 


-2.64 x 10" 1 


32.10 


10- 3 


0.01 


2.92 x 10 1 


1416.2 


333.4 


313.6 


1501.5 


-1.79 x 10 _1 


44.90 


10~ 3 


1.00 


5.47 x 10~ 2 


1.2 


2.2 


0.0 


4.5 


1.67 x 10~ 3 


1.04 


10- 3 


1.00 


6.25 x 10~ 2 


1.4 


3.8 


0.0 


5.5 


2.72 x 10~ 3 


1.06 


10~ 3 


1.00 


9.38 x 10~ 2 


2.1 


11.6 


0.1 


9.8 


8.25 x 10~ 3 


1.14 


10- 3 


1.00 


1.25 x 10 _1 


2.7 


20.3 


0.4 


13.0 


1.39 x 10~ 2 


1.23 


10- 3 


1.00 


1.56 x 10 _1 


3.4 


30.1 


1.1 


16.5 


2.31 x 10~ 2 


1.33 


10~ 3 


1.00 


2.34 x 10 _1 


5.1 


58.3 


4.5 


27.4 


3.81 x 10~ 2 


1.75 


10- 3 


1.00 


3.13 x 10" 1 


6.9 


77.5 


9.4 


36.8 


6.01 x 10" 2 


2.22 


10" 3 


1.00 


4.69 x 10" 1 


10.3 


135.4 


6.6 


78.1 


1.07 x 10" 1 


3.48 


10- 3 


1.00 


6.25 x 10- 1 


13.7 


160.6 


10.0 


107.8 


1.38 x 10" 1 


4.67 


10" 3 


1.00 


7.81 x 10" 1 


17.1 


164.2 


13.2 


137.9 


1.40 x 10" 1 


5.75 


10- 3 


1.00 


9.38 x 10" 1 


20.6 


159.7 


16.4 


165.8 


1.31 x 10" 1 


6.90 


10" 3 


1.00 


1.25 


27.4 


532.9 


61.3 


216.7 


-4.23 x 10" 1 


9.13 


10- 3 


1.00 


1.56 


34.3 


573.4 


73.2 


255.1 


-4.31 x 10- 1 


10.80 


10- 3 


1.00 


3.13 


68.5 


636.4 


114.2 


400.3 


-4.32 x 10" 1 


16.10 


10" 3 


1.00 


4.69 


102.8 


562.3 


140.0 


511.0 


-4.38 x 10" 1 


20.10 


10- 3 


1.00 


6.25 


137.0 


462.6 


140.2 


584.5 


-4.03 x 10- 1 


21.80 


10- 3 


1.00 


9.38 


205.5 


514.7 


214.5 


753.3 


-3.89 x 10" 1 


28.10 


10~ 3 


3.00 


7.92 x 10~ 2 


1.1 


1.3 


0.0 


4.9 


2.17 x 10~ 3 


1.03 


10~ 3 


3.00 


8.80 x 10~ 2 


1.2 


2.8 


0.1 


6.2 


4.57 x 10~ 3 


1.05 


10- 3 


3.00 


1.32 x 10 _1 


1.8 


13.8 


0.3 


12.6 


2.04 x 10~ 2 


1.13 


10~ 3 


3.00 


1.76 x 10 _1 


2.5 


27.8 


0.5 


18.5 


3.65 x 10~ 2 


1.23 


10- 3 


3.00 


2.64 x 10 _1 


3.7 


59.7 


2.1 


36.0 


7.10 x 10~ 2 


1.67 


10~ 3 


3.00 


3.52 x 10 _1 


4.9 


87.3 


3.6 


51.0 


1.04 x 10 _1 


2.06 


10- 3 


3.00 


4.40 x 10 _1 


6.2 


110.5 


5.1 


64.6 


1.24 x 10 _1 


2.42 


10" 3 


3.00 


5.28 x 10" 1 


7.4 


129.8 


6.3 


77.3 


1.45 x 10" 1 


2.74 


10- 3 


3.00 


6.16 x 10- 1 


8.6 


142.8 


7.8 


90.0 


1.59 x 10" 1 


3.07 


10" 3 


3.00 


7.04 x 10" 1 


9.8 


153.8 


8.7 


101.5 


1.65 x 10" 1 


3.32 


10" 3 


3.00 


8.80 x 10" 1 


12.3 


165.9 


10.9 


124.7 


1.66 x 10" 1 


3.87 


10- 3 


3.00 


1.14 


16.0 


459.3 


40.6 


161.6 


-4.40 x 10- 1 


5.44 


10" 3 


3.00 


1.32 


18.5 


501.6 


50.8 


177.7 


-4.62 x 10" 1 


5.98 


10- 3 


3.00 


1.76 


24.6 


565.5 


70.1 


217.6 


-4.99 x 10" 1 


7.19 


10" 3 


3.00 


2.64 


36.9 


633.4 


94.7 


285.1 


-5.18 x 10" 1 


9.00 


10- 3 


3.00 


3.52 


49.2 


669.2 


110.3 


335.3 


-5.25 x 10- 1 


10.30 


10- 3 


3.00 


6.16 


86.1 


682.3 


153.1 


441.0 


-5.39 x 10" 1 


12.20 


10- 3 


3.00 


1.76 x 10 1 


246.1 


753.1 


325.1 


625.8 


-5.41 x 10" 1 


13.80 


10" 3 


5.00 


4.75 x 10" 2 


1.2 


1.5 


0.1 


6.7 


3.79 x 10~ 3 


1.05 


10- 3 


5.00 


7.12 x 10" 2 


1.8 


6.2 


0.6 


11.5 


1.37 x 10" 2 


1.11 


10" 3 


5.00 


9.49 x 10" 2 


2.5 


13.5 


0.9 


18.0 


2.67 x 10" 2 


1.19 


10- 3 


5.00 


1.19 x 10- 1 


3.1 


21.2 


1.6 


27.8 


3.53 x 10" 2 


1.37 


10" 3 


5.00 


1.42 x 10" 1 


3.7 


28.5 


2.4 


34.8 


4.76 x 10" 2 


1.49 


10- 3 


5.00 


1.66 x 10" 1 


4.3 


35.9 


2.9 


40.8 


5.70 x 10" 2 


1.59 


10" 3 


5.00 


1.90 x 10" 1 


4.9 


42.6 


3.1 


46.2 


6.70 x 10" 2 


1.68 



17 



Table A.4: Continued. 
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Table A.4: Continued. 
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